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Scattered  Data  Near-Interpolation  with 
Application  to  Discontinuous  Surfaces 


Renata  Besenghi  and  Giampietro  Allasia 


Abstract.  This  paper  discusses  a particular  type  of  function  approxima- 
tion on  scattered  data  in  a general  number  of  variables,  and  its  application 
to  surface  representation  with  imposed  conditions.  If  the  given  function 
values  are  subject  to  errors,  it  is  not  appropriate  to  interpolate  the  function 
at  the  data  in  the  sense  of  exact  matching.  As  a consequence,  we  formulate 
a weakened  version  of  the  classical  scattered  data  interpolation  problem, 
and  give  a simple  and  efficient  procedure  to  obtain  near-interpolation  for- 
mulas. Near-interpolants  enjoy  many  remarkable  properties,  which  are 
very  useful  from  both  theoretical  and  practical  points  of  view  (shape  pre- 
serving properties,  operator  positivity,  subdivision  techniques,  parallel  and 
multistage  computation).  Applications  of  near-interpolants  to  the  rep- 
resentation of  surfaces,  in  particular  with  faults,  are  discussed  in  detail 
(parameter  values,  localizing  weights,  etc.). 


§1.  Introduction 

In  many  applications,  the  given  function  values  are  subject  to  errors;  hence 
it  is  not  appropriate  to  interpolate  the  function  at  the  data  in  the  sense  of 
exact  matching,  but  it  seems  more  appropriate  to  approximate  the  function 
or,  more  precisely,  to  get  a relaxed  interpolation  or  near-interpolation.  Data 
requiring  near-interpolation  by  scattered  data  methods  occur  in  virtually  ev- 
ery field  of  science  and  engineering.  Sources  include  both  experimental  results 
(experiments  in  chemistry,  physics,  engineering)  and  measured  values  of  phys- 
ical quantities  (meteorology,  oceanography,  optics,  geodetics,  mining,  geology, 
geography,  cartography),  as  well  as  computational  values  (e.g.,  output  from 
finite  element  solutions  of  partial  differential  equations). 

As  a consequence  of  this  remark,  we  formulate  a relaxed  version  of  the 
classical  multivariate  interpolation  problem  at  scattered  data  points. 
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Definition  1.  Given  a set  of  points  Sn  = { Xi , i = distinct  and 

generally  scattered,  in  a domain  D C IRs,(s  > 1),  with  associated  values 
{fi,  i = 1, . . . ,n},  and  a linear  space  <fr(D)  spanned  by  continuous  real  basis 
functions  gj(x;r ) with  x € D,  r > 0,  and  j = the  multivariate 

near-interpolation  problem  at  scattered  data  consists  in  finding  a function 
F{x\r)  € such  that 

n 

F(xiir)-^2ajgj(xi;r)  = fi  + ei(r),  i = l,...,n,  (1) 

i= i 

and 

limei(r)  = 0.  (2) 

We  observe  that  r works  as  a parameter,  and  the  limit  case  of  F{x]  r)  when  r 
vanishes 

F(x)  = F(a;;0)  = lim  F(x;r) 

is  an  interpolation  operator.  If  F{x\r)  is  specified,  then  the  Ci{r)  in  (1)  are 
known;  these  near-interpolation  errors  at  the  nodes  must  not  be  confused 
with  the  unknown  errors  which  affect  the  corresponding  function  values  /,. 
However,  it  is  reasonable  to  get  things  so  that  the  e,(r)  and  the  errors  on  /; 
are  quantities  of  the  same  order. 

In  Section  2 we  give  a constructive  procedure  to  obtain  a wide  class  of 
near-interpolation  formulas.  These  enjoy  many  interesting  properties  which 
are  listed  in  Section  3.  A crucial  point  in  near-interpolation  is  the  proper 
choice  of  the  parameter  r in  (1)  and,  eventually,  of  other  parameters;  the 
matter  is  discussed  in  Section  4.  Finally,  Section  5 is  devoted  to  the  application 
of  near-interpolation  to  modelling  faults. 

§2.  Construction  of  Near-Interpolants 

To  solve  the  classical  interpolation  problem,  one  can  consider  basis  functions 
which  depend  on  the  nodes  and,  moreover,  are  cardinal.  The  method  of  car- 
dinal basis  functions  involves  selecting  continuous  cardinal  functions  gj  : D — > 
IR,  (j  = 1 ,...,n),  such  that  gj{  Xi)  = 6ij  , (i  = 1,  ...,n),  where  Sij  is  the 
Kronecker  delta  operator,  and  setting  up  the  interpolation  operator  F in  the 
form 

n 

Fix ) = ][]  fi  9j{x)- 
i= i 

The  corresponding  near-interpolation  problem  considers  basis  functions 
gj(x\r),  which  are  no  longer  cardinal,  but  gj(x-,r)  — > Qj(x)  for  r — > 0.  If  such 
gj(x\r)  are  given,  then 

n n 

F{x;  r)  = fi  9iix 5 r)  = F{x)  + ]T  fj  [ gj{x ; r)  - ^(x)]  (3) 

i= i i= i 

represents  a solution  of  the  near-interpolation  problem.  In  this  relation  the 
terms  e,(r)  = F(x,;  r)  — fi,  (i  = 1, . . . , n),  are  uniquely  determined  and  satisfy 
(2). 
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As  an  example,  let  us  consider  the  near-interpolant  Shepard’s  formula  in 
the  product  form 


, A lTk=i,kfr[d2(x’xi<)  + rf 

x’ r)  U ELi  UU **) + r]t>  ’ 

where  d(x,y)  is  the  Euclidean  distance  between  x and  y.  and  (5  > 0,  or  in  the 
equivalent  barycentric  form 


n 

p(x;  r)  = £ fi 
j= 1 


[d2(x,Xj)  + r]  & 

Y,Z=M'2(x,xh)  + r}-P  ’ 


p(xi-0)  = /*, 


i = 1, . . . , n.  (5) 


This  formula  no  longer  interpolates  for  r > 0,  but  p(x\  r)  — > p(x)  as  r — * 0, 
where  p(x)  is  the  well-known  Shepard’s  formula  [8,  1]. 

Examining  the  structure  and  the  basic  idea  of  Shepard’s  operator  suggests 
a simple  and  efficient  procedure  to  obtain  an  interpolation  formula  [1].  The 
corresponding  way  to  obtain  a near-interpolation  formula  is  contained  in 

Definition  2.  Let  a(x,y,r),  with  x,y  € D and  r > 0,  be  a continuous 
positive  real  function  such  that 


lim  a(x,y;r)  = a(x,y), 


(6) 


where  a(x,y)  > 0 , if  x ^ y and  a(x,y)  = 0 , if  x = y for  all  x,y  S D.  Define 
now  the  functions  gj(x;r)  by  the  equations 

II  <*(*.**;*•) 

Eh=i  nfc=i,ife#h  «(*.**;*■) 

and  the  near-interpolant  F(x;  r)  by 


E(x;  r)  = £ /,  *(«; r)  = £ /,  . 

j-l  j=i  2-fh=l  llk=l,k^ha\x’xk’ 


r ) ’ 


(8) 


or  equivalently  by 


F(X;r)  ~ ELit /aix’Zr)  ’ F{Xi’0)  ~ fi’ 


i = 1, . . . ,n.  (9) 


Many  choices  are  possible  for  the  function  a(x,y;r)  in  (6);  there  are  no 
constraints  engendered  by  the  set  Sn,  that  is,  the  distribution  of  the  nodes 
is  irrelevant.  Nevertheless,  experience  suggests  identifying  a with  a radial 
function 


a(x,y;r)  = cj)(\\x  - y\\2  + r), 
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where  ||  - 1|  is  a convenient  norm.  As  an  example,  choosing  the  Euclidean  norm 
||  ■ ||2  and 

a(x,y;r)  = (\\x-y\\l+rf,  p > 0,  (10) 

we  obtain  from  (8)  the  near-interpolant  Shepard’s  formula  (4). 

It  is  often  convenient  to  consider  the  function  1 p(t-  r)  : R>o  — > It,  as- 
sociated with  (f>(\\x  — y\\2  + r)  and  defined  as  <p(£;r)  = (f>(t2  + r).  We  have 
just  seen  </?(£;  r)  = (t2  + r)@  in  (10);  another  possible  expression  that  works  is 
<p(f;  r)  = (f2  + r )&  exp(7  t2),  (7  > 0).  Considering  several  functions  <p  can  be 
useful  in  order  to  compare  their  behaviour  and  choose  the  most  suitable  for 
use. 


§3.  Properties  of  Near-Interpolants 

Near-interpolants,  as  given  in  Definition  2,  enjoy  many  interesting  properties. 
We  list  some  of  them. 

A)  The  near-interpolant  F(x;r)  in  (8)  or  (9)  is  a weighted  arithmetic  mean 
of  the  values  /j,  ( j = 1 ,...,n),  since  0 < gj(x;r)  < 1 and  Y^j=\9j(x>r)  — 
1.  As  a consequence  F(x;r)  satisfies  the  betweeness  property  min;/;  < 
F(x;  r)  < max;  /; , and  reproduces  exactly  any  constant  function  f(x)  = c , 
that  is,  if  /;  = c,  (i  = 1, . . . ,n),  then  F(x\r)  = c.  Moreover,  F(x;r), 
considered  as  a functional  on  the  set  of  functions  / : D — > 1R,  is  linear  and 
positive. 

B)  If  a(x,y;r)  is  infinitely  differentiable  with  respect  to  the  pth  compo- 
nent of  x = (xM, . . . , x^)  for  all  x,y  G D and  p = 1, . . . , s,  then  F(x;  r)  is 
also  infinitely  differentiable  with  respect  to  the  xl-p'> . For  example,  choosing 
a(x,y;r)  as  in  (10),  F(x;r ) = p(x\r)  in  (4)  can  be  differentiated  as  many 
times  as  desired. 

C)  If  a is  a radial  function,  F(x;  r)  enjoys  some  properties  of  invariance  with 
respect  to  affine  transformations.  In  particular,  with  the  Euclidean  norm,  we 
have  that  F(x\  r)  is  invariant  under  translation  and  rotation,  but  not  scalar 
invariant. 

D)  Subdivision  techniques  can  be  applied  to  near-interpolants  achieving  re- 
markable results,  very  well  suited  for  parallel  computation  [3].  Let  us  make  a 
partition  of  the  set  Sn  on  the  domain  D into  q subsets  Snj , so  that  the  jth 
subset,  (j  = 1, . . . , q),  consists  of  the  nodes  Xj  1 , Xj2,  ■ ■ ■ > xjnj , with  n\  + 712  + 
• • • + nq  = n,  and  the  values  fjkj  ,{j  = 1, . . . , q\  kj  = 1, . . . , rij),  correspond  to 
the  nodes  Xjkj  ■ The  indexing  of  the  nodes  in  the  subsets  may  not  depend  on 
the  indexing  in  the  set,  provided  the  biunivocity  is  saved. 

Given  Snj  = {xji,  xj2, . . . , xjnj },  (j  = 1,  — , 9),  let  Sn  = Sni  U5„2  U-  • -U 
Snp  and  Snq  fl  5Ur  = 0 for  q ^ r,  then  F(x;  r)  in  (9)  can  be  rewritten  in  the 
form 


9 A • nj 

Fsn(x\r ) = Y^FSnj{x\r)  q 3 , where  Aj  = ^ l/a{x,xjkj;r).  (11) 

j=i  L,j=iAi  fcj=i 
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E)  As  a consequence  of  (11)  the  following  multistage  procedure  works  very 
well.  In  the  first  stage,  a given  set  of  nodes  Sni  = {xi,i  = 1, ...,ni}  is 
considered,  and  the  corresponding  near-interpolant  Fsni(x;r)  is  evaluated. 
In  the  second  stage,  it  is  required  to  enlarge  the  considered  set  Sni,  taking 
the  union  of  it  and  another  set  of  nodes  5„2  = {xj,j  = 1, . . . ,712}.  Now  the 
near-interpolating  function  referred  to  the  union  set,  i.e.,  Tsniusn  (x< r)i  with 
Sni  l~l  Sn2  = 0,  can  be  obtained  simply  by  evaluating  the  near-interpolant 
Fs„2  (x;  r),  corresponding  to  the  added  set  S»2,  and  using  the  relation 


FSniuSn2  (x'r)  = 


Fsni  (x;  r)A1  + Fs„2  (a;;  r)A2 

Ai  + A2 


(12) 


where  Ai  = ^la{xixi\r)^  -^2  = l/a(x>xjir)t  and  A\  is  known. 

The  procedure  can  be  repeated  as  many  times  as  required. 


F)  Near-interpolants  have  the  remarkable  property  that  an  additional  node, 
say  xn+\,  can  be  added  to  the  interpolation  set  Sn  by  simply  combining  an 
extra  term  with  the  original  formula.  The  goal  is  achieved  by  using  a particular 
case  of  (12),  that  is,  the  recurrence  relation 

r,  , \ FSn(x-,r)An  + fn+1  l/a{x,xn+v,r) 

= A.  + !/«(»,  ' 

where  = X)fc=i  1 /a(x,xk;r). 

G)  It  is  often  convenient  to  extend  (8),  or  better  (9),  in  the  following  way: 


n 


Fi{x\r)  = J2fi 

3- 1 


r(x,xj-,j)  l/a(x,Xj;r) 
T,h=iT(x’xh;i)  l/a(x,xh;r)  ’ 


(13) 


where  t(x,  y\ 7),  with  x,y  € D and  7 > 0,  is  a continuous  positive  real 
function.  Choosing  suitably  r(x,y,  7),  one  can  modify  the  weights  in  (13) 
in  order  either  to  cancel  a useless  characteristic,  or  to  introduce  a new  fea- 
ture. In  particular,  it  is  possible  to  localize  the  method  considering  a factor 
r(a;,i/;7)  rapidly  decreasing  with  distance  [2],  The  formulas  obtained  in  this 
way  maintain,  in  general,  the  analytical  and  computational  properties  of  the 
corresponding  original  ones. 

The  use  of  the  exponential-type  function 

x{x,y\l)  =exp(— ■ -)\\x-y\\l)  (14) 


is  suggested  by  McLain  [4]  for  Shepard’s  formula;  he  observes  that  much  more 
accurate  results  can  be  obtained  in  this  way.  The  use  of  exponential-type 
weights  increases  the  computational  effort,  but  generally  this  drawback  can 
be  tolerated. 

The  value  of  the  parameter  in  the  mollifying  function  r( x,  y\ 7)  may  de- 
pend on  the  nodes,  as  happens  in  the  popular  case  [4] 

T(x,Xj-,pj)  = f 1-— — , (15) 

where  pj  is  the  radius  of  the  circle  of  support  at  the  point  Xj,  and  (u)+  > 0 
if  u > 0,  (u)+  = 0 if  u < 0. 
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H)  The  precision  of  the  operator  F(x;r)  can  be  increased,  considering  the 
Taylor  expansion  for  the  function  / in  each  node  Xj  instead  of  the  function 
value  fj.  This  leads  to  the  following  extension  of  (9). 

If  / € Cm(D)  and  Tj(x ) is  the  truncated  Taylor  expansion  for  / up  to 
derivatives  of  order  m evaluated  at  the  point  Xj  and  referred  to  the  displace- 
ment hj  = x — Xj,  with  hj  C D,  then  the  operator 


n 

F2{x-,r)  = Y^Tj(x)  9j(x\r),  (16) 

3 = 1 

near- interpolates  to  Tj{x)  at  x = Xj.  In  this  form,  1*2(2:;  r)  reproduces  exactly 
algebraic  polynomials  of  degree  < m. 

Combining  the  modifications  in  (13)  and  (16),  we  have 


F3(x-,r)  = Y/Tj{  x) 
3= 1 


T(x,Xj;-y)  l/a(x,Xj-,r) 
£h= iT(x,xh]l)  1 /a{x,xh;r)  ' 


(17) 


Obviously,  the  technique  calls  for  additional  derivative  values  that  are  not 
normally  available  as  data.  A more  practical  solution  is  discussed  below. 

K)  For  simplicity,  we  refer  here  to  an  Euclidean  radial  function  a(x,  y,  r ) = 
<j>(\\x  — y Hi  + r),  because  in  this  case  the  procedure  is  well  established.  The 
primary  modifications  required  involve  using  t(x,  y;  7)  to  localize  the  overall 
approximation,  and  replacing  fj  with  a suitable  “local  approximation”  to 
the  surface.  To  carry  out  the  approximation  (17),  a practical  way  is  to  get, 
in  a first  stage,  local  approximants  Mj(x)  to  f(x)  at  the  points  Xj,(j  = 
1 ,...,n),  obtained  by  means  of  the  moving  weighted  least-squares  method 
using  weight  functions  with  reduced  compact  support.  Then,  in  a second 
stage,  the  near-interpolating  operator  is  expressed  as  a convex  combination  of 
the  local  approximants 


F4(x;r)  = Mj(x) 
i- 1 


r{x,Xj\i)  l/<f>{\\x-Xj\\l  + r) 

£h=!  T(x’xh-,l)  IMII®  - ajftlli  +r)  ' 


(18) 


In  particular,  by  (18),  (10),  and  (14), 


n 

Mi  (*;r)  = YMAX) 

3= 1 


exp(— 7||x  — 1 1 2 ) i\\x  - Xj\\l  + r)  £ 

ELi  exp(-7||x  - 27,111)  (II*  ~ */>lll  + r)~p  ’ 


(19) 


which  extends  (5). 

Very  good  performance  is  achieved  by  a version  of  (18)  which  uses  quad- 
ratic approximations  for  Mj(x),  and  mollifying  functions  given  by  (15).  This 
method  has  been  developed  by  Franke  and  Nielson  [4],  and  Renka  [7]  for 
Shepard’s  operator. 
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§4.  Determining  Parameter  Values 

The  near-interpolating  operator  (19)  works  very  well  in  a large  variety  of 
cases.  Our  attention  is  here  focused  on  finding  values  of  the  parameters  r,  (3 
and  7,  which  can  be  regarded  as  “optimal”  from  a practical  viewpoint.  The 
considerations  which  follow  are  mainly  based  on  experiments. 

A relatively  small  increase  of  the  parameter  r in  (19)  in  some  right  neigh- 
bourhood of  zero  has  a considerable  effect  on  the  behaviour  of  pi(x\ r).  In  fact, 
if  r is  small,  x fixed  and  near  to  the  node  Xj* , then  the  value  of  pj*  (x;  r)  equals 
nearly  one;  but,  if  r increases,  gj-(x;  r)  decreases.  Since  )T™=1  9j{x\ r)  = 1> 
diminishing  of  gj-(x;r)  makes  the  other  weight  values  gj(x;r),  j yt  j*f  in- 
crease. Summing  up,  if  the  weight  attributed  to  /j*  in  pi(a :;r)  decreases, 
then  pi{xj+ ; r)  = /j.  -f  (r)  diverges  from  fj .,  namely  tj- (r)  increases  and 
reduces  the  accuracy  of  pi  (xj*\r). 

Introducing  the  parameter  r in  (19),  and  in  particular  in  (5),  has  the 
effect  that,  in  general,  the  gradient  of  the  rendered  surface  is  not  zero  at  the 
nodes.  As  a consequence,  the  surface  is  considerably  smoother  than  for  r — 0. 
However,  if  r is  too  small,  the  first  derivatives  of  pi(x',  r)  are  highly  oscillating 
and  their  values  are  nearly  zero.  Clearly,  the  goal  is  to  choose  an  “optimal” 
value  of  r,  such  that  pi(x;  r)  does  not  exhibit  the  characteristic  irregularities 
of  the  basic  Shepard’s  formula,  but  at  the  same  time,  it  maintains  a sufficient 
computational  accuracy,  in  particular  at  the  nodes. 

The  search  for  the  optimal  value  of  r can  be  done  by  many  applications 
of  (19)  with  different  values  of  the  parameter,  and  then  by  choosing  that  value 
which  minimizes  the  global  root  mean  square  error.  Although  this  is  currently 
considered  in  the  literature,  the  estimate  of  r is  not  a simple  matter;  in  a 
sense,  it  can  be  compared  with  the  analogous  difficult  problem  of  computing 
the  optimal  value  of  the  parameter  in  multiquadric  interpolants. 

The  optimal  value  of  the  parameter  (3  has  been  determined  with  particular 
attention  to  computational  accuracy.  The  performance  analysis  on  some  test 
functions  proposed  by  Franke  leads  to  prefer  the  value  (3  = 3/2. 

As  for  the  optimal  value  of  the  parameter  7 in  the  strongly  localizing 
function  (14),  McLain  has  proposed  7 = 1.62n/diam(£>),  where  n is  the  num- 
ber of  nodes  and  diam(D)  is  the  diameter  of  D.  However,  this  value  is,  in 
general,  too  large,  whereas  it  is  sufficient  to  consider  for  7 a value  of  the  order 
of  tens. 


§5.  Application  to  Modelling  Faults 

Using  the  near-interpolating  operator  pi(x;  r)  of  (19),  with  a suitable  value  of 
the  parameter  r,  instead  of  the  corresponding  interpolating  operator  pi(x;  0), 
increases  considerably  the  performance  of  the  approximation  in  a rich  variety 
of  applications,  because  it  permits  consideration  of  supplementary  information 
connected  with  the  characteristics  of  the  examined  problem.  A typical  case 
occurs  with  surface  discontinuities,  in  particular  faults,  which  are  frequently 
met  when  modelling  geological  surfaces. 
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Following  Shepard  [8],  we  observe  that,  if  some  physical  barrier  such  as  a 
fault  separates  the  set  of  nodes,  the  relationship  between  nodes  on  the  opposite 
sides  of  the  barrier  may  be  attenuated.  Through  the  inclusion  of  barriers,  a 
user  may  specify  discontinuities  in  the  metric  space  in  which  the  distance 
between  two  points  is  calculated  to  simulate  this  attenuation.  Suppose  a 
“detour”  of  length  b(x,Xk)  were  required  to  go  over  the  barrier  between  x, 
the  current  near-interpolation  point,  and  the  node  x k-  The  quantity  b(x,Xk) 
is  considered  the  strength  of  the  barrier,  and  an  effective  distance  between  x 
and  Xk  is  given  by 


d*(x,xk)  = y/[d(x,Xk)]2  + \b(x,xk)}2. 

This  definition  is  general  so  that  if  no  barrier  separates  x and  aq , then 
b(x,  Xk)  = 0 and  d*(x,  Xk)  — d(x,  Xk )•  Because  of  the  discontinuity  in  effective 
distance  as  the  near-interpolation  point  x crosses  the  barrier,  the  rendered 
surface  will  be  discontinuous  at  the  barrier. 

Since  extensive  tests  [4]  have  shown  that  the  modified  quadratic  Shepard’s 
method  performs  very  well  for  a variety  of  data  sets,  Franke  and  Nielson  [5] 
have  chosen  it  as  a basis  to  investigate  the  problem  of  simulating  faults.  Our 
approach  uses  instead  the  near-interpolant  (19),  with  significant  differences  in 
distance  penalty,  localizing  functions,  fault  forms,  etc.,  as  compared  to  Franke 
and  Nielson. 

The  possibility  of  having  to  model  faults  can  occur  in  different  ways  [5]; 
to  save  space,  we  limit  our  attention  to  the  following  case:  there  is  a known 
fault  line  T C D C R2,  in  a known  location,  with  a known  jump.  More 
complicated  situations  (see,  e.g.,  [5,6])  require  extensive  considerations  that 
will  be  discussed  in  a further  work. 

As  a first  step,  it  is  convenient  to  focus  on  the  basic  situation  in  which 
the  fault  line  T is  a segment  l and,  moreover,  the  jump  is  constant  along  l. 
Then,  a known  polygonal  curve  can  be  considered  as  a fault  line;  in  fact,  the 
reduction  to  the  case  of  a fault  line  segment  is  straightforward  by  using  the 
subdivision  procedure  considered  in  Section  3.  In  principle,  any  curve  can  be 
considered  as  a fault  line,  provided  it  is  well  approximated  by  a polygonal. 
Another  extension  consists  in  considering  a jump  varying  along  the  fault  line. 
Also  the  reduction  to  the  basic  case  is  now  possible,  subdividing  the  fault  line 
into  a convenient  number  of  segments,  and  using  a mean  value  of  the  jump 
for  each  segment. 

To  deal  with  the  basic  case,  we  modify  the  value  of  the  parameter  r in 
(19)  in  order  to  take  the  jump  into  account.  Let  x be  the  near-interpolation 
point,  Xk  a node  and  l*  the  segment  joining  x and  Xk-  Then  for  x,Xk  I we 
set 

f ropt , if  i n r = 0, 

l b(x,xk),  if  l n /*  ^ 0, 

where  the  quantity  ropt  is  the  optimal  value  obtained  for  r in  (19)  on  the 
opposite  sides  of  the  fault  and  b(x,Xk)  represents  the  “effort”  required  to  go 
over  the  barrier,  due  to  the  discontinuity  dividing  the  two  points.  If  the  jump 


Fig.  1.  Function  fi(x,  y):  near-interpolation  and  signed  error  surfaces. 


Fig.  2.  Function  f2(x,y)\  near-interpolation  and  signed  error  surfaces. 


is  constant  or  almost  constant  along  l,  it  is  possible  to  simply  set  b(x,  Xf.)  — h, 
where  h is  the  jump  size. 

Formula  (19),  after  these  adjustments  in  the  parameter  r,  gives  results 
quite  good  both  for  the  appearance  of  the  graphic  representation  and  the  accu- 
racy in  computation.  Comparing  the  rendered  surface  with  the  one  obtained 
by  the  modified  quadratic  Shepard’s  formula  shows  that  the  introduction  of 
the  parameter  r gives  a smoother  surface  which  is  closer  to  the  approximated 
function. 

Our  procedure  has  been  used  to  fit  the  test  function  proposed  by  Franke 
and  Nielson  [5]  using  their  set  of  nodes.  Numerous  tests  were  also  made  on 
other  surfaces.  We  present  two  examples  of  the  rendered  surfaces  and  the 
signed  error  surfaces  for  the  functions 


if  0 < x < 0.5, 
if  x > 0.5; 


h{x,y) 


y - x + 1,  if  0 < x < 0.5, 

0,  if  0.5  < x < 0.6, 

0.3,  if  x > 0.6, 
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defined  on  the  unit  square  (see  Fig.  1 and  Fig.  2).  We  used  the  parameter 
values  r = 0.0036,/?  = 1.5,7  = 24,  and  r = 0.0025,/?  = 1.5,7  = 30  respec- 
tively, and  once  again  the  set  of  nodes  of  Franke  and  Nielson.  The  errors  can 
be  considerably  reduced  by  adding  more  information  on  the  faults;  in  fact, 
the  employed  set  of  nodes  is  not  obviously  an  ad  hoc  choice. 
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